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ABSTRACT 

This paper describes a finite volume based numerical algorithm that allows multi-dimensional 
computation of fluid flow within a system level network flow analysis. There are several 
thermo-fluid engineering problems where higher fidelity solutions are needed that are not within 
the capacity of system level codes. The proposed algorithm will allow NASA’s Generalized 
Fluid System Simulation Program (GFSSP) to perform multi-dimensional flow calculation 
within the framework of GFSSP’s typical system level flow network consisting of fluid nodes 
and branches. The paper presents several classical two-dimensional fluid dynamics problems 
that have been solved by GFSSP’s multi-dimensional flow solver. The numerical solutions are 
compared with the analytical and benchmark solution of Poiseulle, Couette and flow in a driven 
cavity. 


INTRODUCTION 

In thermo-fluid engineering applications, system level codes are typically used to find flow and 
pressure distribution in a complex flow network. On the other hand, Navier-Stokes codes are 
used when detailed knowledge about the flow is needed for design or investigate a failure 
scenario. System level codes are often run independently to provide boundary conditions for 
Navier Stokes codes. There has not been much success in the integration of these codes to 
perfonn any coupled analysis. However, there are situations where integrated analysis brings 
value to the design. One such example is the propellant feed to a rocket engine from a stratified 
cryogenic tank. While the bulk of the feed system analysis can be perfonned by a system level 
flow network code, the stratification is a multi-dimensional phenomenon and requires Navier- 
Stokes analysis in conjunction with network flow analysis which models pressurization and 
venting. 
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In order to analyze such problems, it was decided to expand NASA’s Generalized Fluid System 
Simulation Program (GFSSP) [1] to account for multidimensional flow effects. GFSSP is a 
finite volume based [2] network flow analysis code where a complex flow network is discretized 
into nodes and branches. The mass, energy and species conservation equations are solved at the 
nodes while momentum equations are solved at the branches. The one dimensional momentum 
equation in a flow network typically consists of pressure gradient, friction and fluid acceleration 
tenns. In this paper the momentum equation has been extended to include shear due to velocity 
gradient and transport of longitudinal momentum due to transverse velocity component from 
neighboring branches in order to account for multi-dimensional effect of the flow. The 
unstructured nature of the code has been preserved so as to retain flexibility for irregular 
geometries. 

MATHEMATICAL FORMULATION & VERIFICATION OF RESULTS 


GFSSP employs an upwind differencing scheme to derive the conservation equations of mass, 
momentum, energy and species conservation and solves the system of equations by a 
combination of successive substitution scheme and simultaneous (Newton-Raphson) scheme [3], 
This section describes the extension of one-dimensional momentum equation to include a) the 
shear term in laminar flows due to velocity gradient, and b) transport of longitudinal momentum 
due to transverse velocity component from neighboring branches. 


Shear term in momentum equation 


The two-dimensional, steady-state Navier-Stokes Equations in rectangular Cartesian co-ordinate 
can be written as: 
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The first tenn on the left hand side of equation (1) corresponds to the longitudinal acceleration 
(inertia) term within GFSSP. The second tenn on the left hand side of equation (1) corresponds 
to the transverse momentum exchange. The first tenn on the right hand side of equation (1) 
conesponds to the pressure tenn within GFSSP. The second tenn on the right hand side of 
equation (1) corresponds to the gravity term within GFSSP. The third tenn on the right hand 
side of equation (1) is negligible (based on an order of magnitude argument). The fourth tenn on 
the right hand side of equation (1) corresponds to the friction tenn of GFSSP. It is this fourth 
tenn that represents the shear term that is sought. 

The discretized form of shear tenn is expressed in equation 2. 
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Figure 1, below, represents a set of nodes and branches for which shear forces are exchanged. 
Let branch 12 represent the branch for which the shear force is to be calculated. Branches N12 
and S12 represent the parallel branches which will be used to calculate the shear force on branch 
12. Let YS be the distance between branches 12 and S12, and let YN be the distance between 
branches 12 and N12. Let AS be the shearing area between branches 12 and S12, while AN is 
the shearing area between branches 12 and N 12. 
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Figure 1 : Grid for Shear Exchange 


It is possible (and highly probable) that YS & YN will not be the same nor exactly parallel to 
branch 12, therefore, a modified differencing scheme must be used. Equation (3) below is the 
proposed differencing scheme to be implemented into GFSSP: 
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Equation (3) can be generalized to n-number of parallel branches at any position around branch 
12 as shown in equation (4). 
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where the summation subscript i represents the i th parallel branch to branch 12. 

The expression for the shear effect of the wall on branch 12 is given in equation (5). 
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If there were multiple walls adjacent to branch 12, then equation (5) could be generalized into 
equation (6). 
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Equation (7) represents the laminar shear fonnulation where i is the current branch, np; is the 
number of parallel branches to branch i, and nsj is the number of parallel solid walls to branch i. 
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Let’s examine two situations which will benchmark the laminar shear term of equation (7). The 
first situation is that of Poiseulle flow between infinite flat plates. The second situation is that of 
Couette flow between infinite parallel flat plates. The implementation of equation 7 in GFSSP 
has been verified by comparing its solution with two bench mark laminar flow examples: 
Poiseulle Flow and Couette Flow 

POISEULLE FLOW 


Consider the flow between two fixed flat plates shown below in Figure 2. The flow is pressure 
driven and assumed to be fully developed. The length of the plates is 1000 inches, they are 1 
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inch apart, the fluid density is 12 lb m /ft , the viscosity is 1 lb m /(ft sec), and the pressures at the 
upstream and downstream ends are 20 psi and 10 psi, respectively. Figure 3 shows an 
approximate velocity profile for this situation. 
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Figure 2: Poiseulle Flow Physical Situation 
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Figure 3: Poiseulle Flow Velocity Profile 


The analytical solution of velocity for Poiseulle flow [4] for situation shown in figure 2 can be 
expressed as: 


u = 0.005(y-y 2 ). 


( 8 ) 


GFSSP model of Poiseulle Flow 

A simple 3 node, 10 branch model (two sets of 5 parallel branches) was constructed to model the 
physical situation described above. The model is shown schematically in Figure 4 below. 
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Figure 4: GFSSP Model of Poiseulle Flow 

Each of the 10 branches (121-125, 231-235) are 500 inches in length. It is assumed that the 
shear area for each branch-branch & branch-wall was 500 in and that the shear distance was 0.2 
inches between adjacent branches and 0.1 inch between walls and their adjacent branches. The 
pressures at nodes 1 and 3 are 20 psi and 10 psi, respectively. A constant density option was 
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used and the corresponding density and viscosity used were 121b m /ft and 1.0 lb m /(ft'sec), 
respectively. Resistance option -02 was used in the initial flow field calculation (for the first 
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Newton-Raphson iteration, after which the shear will replace the friction factor calculation) for 
each of the branches. The bottom and top walls are fixed. 

Figure 5 below shows a comparison between the velocity profiles for the analytical solution and 
the GFSSP 3 node, 10 branch (5 parallel branch) model. As can be seen in Figure 5, the results 
of this crude GFSSP model compare very favorably with the analytical solution. 



Figure 5: Poiseulle Flow Velocity Distribution 


COUETTE FLOW 


Consider the flow between two fixed flat plates shown in Figure 6. The flow is shear driven and 
assumed to be fully developed. The length of the plates is 1000 inches, they are 1 inch apart, the 
fluid density is 12 lb m /ft , the viscosity is 1 lb m /(ft'sec). The pressures at the upstream and 
downstream ends are both set at 10 psi. The bottom plate is fixed, and the top plate is moving 
with a velocity of 100 ft/sec. Figure 7 shows an approximate velocity profile for this situation. 
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Figure 6: Couette Flow Physical Situation 



Figure 7: Couette Flow Velocity Profile 

The analytical solution for Couette flow situation [4] shown in Figure 6 can be expressed as 

u = lOOy . (9) 


GFSSP Model of Couette Flow 

A simple 3 node, 10 branch model (two sets of 5 parallel branches) was constructed to model the 
physical situation described above. The model is shown schematically in Figure 8 below. 



Figure 8: GFSSP Couette Flow Model 

Each of the 10 branches (121-125, 231-235) is 500 inches in length. It is assumed that the shear 
area for each branch-branch & branch-wall was 500.0 in and that the shear distance was 0.2 
inches between adjacent branches and 0.1 inch between walls and their adjacent branches. The 
pressures at nodes 1 and 3 are both 10.0 psi. A constant density option was used and the 
corresponding density and viscosity used were 12.0 lb m /ft and 1.0 lb m /(ft'sec). Resistance option 
-02 was used in the initial flow field calculation (for the first Newton-Raphson iteration, after 
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which the shear will replace the friction factor calculation) for each of the branches. The bottom 
walls are fixed, and the top walls are moving at a velocity of 100 ft/sec. 

Figure 9 below shows a comparison between the velocity profiles for the analytical solution and 
the GFSSP 3 node, 10 branch (5 parallel branch) model. As can be seen in Figure 9, the results 
of the GFSSP model compare nearly identically with the analytical solution. 



Figure 9: Couette Flow Velocity Distribution 

Transverse Momentum Term in Momentum Equation 

The second term on the left hand side of equation (1) corresponds to the transverse momentum 
exchange. The discretized fonn of this term is expressed in equation 10. 

fipvu (pvA u) 

~dy~ V ~ -^^( Ax X A y)( Az ) = (A x)(A z)(pvAu) = m trans Au (10) 

Figure 10, below, represents a set of nodes and branches for which transverse momentum 
exchange will take place. Let the branch 12 represent the current branch which will receive 
transverse momentum from the surrounding branches. Branch S12 represents an adjacent 
parallel branch, while branches SI and S2 represent the adjacent nonnal branches. 
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Figure 10: Grid for Transverse Momentum Exchange 

Now, we will examine the fonnulation for calculating the transverse momentum term for branch 
12. First, calculate the average mass flow rate for the adjacent normal branches: 
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Now, examining Figure 10 we see that a positive transverse mass flow rate is defined as flow 
into the nodes corresponding to the branch in question. Based on this definition of transverse 
mass flow rate, let’s calculate the transverse momentum term: 
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The parallel branch (S12) will contribute to the transverse momentum term to branch 12 when 
m 12 > 0 (since a positive transverse flow rate begins at S12 and ends at 12), and should have a 
negligible contribution when m 12 < o (he. transverse flow rate begins at 12 and ends at S12). 
Equation (13), below, is the upwinding representation of the transverse momentum tenn for 
branch 12. 
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Equation (13) can be generalized for m-parallel branches around branch i, each with an angle Q, 
with respect to branch i, and n^ corresponding transverse connecting branches, each transverse 
branch with an angle of Qjk with respect to branch i. Equation (14) represents this generalized 
version of equation (13). 
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Equation (14) represents the transverse momentum formulation where i is the current branch for 
which transverse momentum is being calculated, m; is the number of parallel branches which 
will be used to calculate transverse momentum, and n;j is the number of connecting transverse 
branches between the current branch (i), and the j th parallel branch. 

The implementation of Equation (14) in GFSSP has been verified by comparing its solution with 
benchmark solution of Burggraf for driven cavity [5]. 

Driven Cavity Flow 

In a square cavity the flow is induced by shear interaction at the top wall as shown in Figure 11. 
The length of each wall is 12 inches. The density of the fluid is assumed constant at 1.00 lb m /ft , 
and the viscosity of the fluid is assumed to be 1.00 lb m /(ft sec). The bottom and side walls are 
fixed. The top wall is moving to the right at constant speed of 100 ft/sec. The corresponding 
Reynolds number for this situation is Re = 100. 



T 


12 inches 


Figure 1 1 : Flow in a Shear Driven Square Cavity 

GFSSP Model of Driven Cavity 

The GFSSP model of the driven cavity consists of 49 nodes (48 of which are internal) and 84 
branches. For numerical stability, one node (node 1) was assigned to be a boundary node with 
an arbitrary pressure of 100 psi. A unit depth (1 inch) was assumed for the required areas. The 
model is shown schematically in Figure 12. 
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It is assumed that the shear area for each branch is 2.0 in". The shear distance between adjacent 
branches is 1.71429 inches and the shear distance between walls and their adjacent branches is 
0.85714 inches. For transverse momentum, only adjacent parallel branches were considered, and 
only connecting branches not associated with the boundary node were used. A constant density 
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option was used and the corresponding density and viscosity used were 1.0 lb m /ft and 1.0 
lb m /(ft'sec), respectively. Resistance option -02 was used in the initial flow field calculation (for 
the first Newton-Raphson iteration, after which the shear will replace the friction factor 
calculation) for all of the branches. The bottom and side walls are fixed. The top walls are 
moving to the right at 100 ft/sec. All parallel angles are 0°, and all transverse angles are 90°. 

Figure 13 shows a comparison between the benchmark numerical solution and GFSSP 7x7 node 
model velocity profiles along a vertical plane at the horizontal midpoint. As can be seen in 
Figurel3, the results of this crude GFSSP model compare very favorably with the benchmark 
numerical solution of Burggraf [6]. 


Top Wall Top Wall Top Wall Top Wall Top Wall Top Wall 



Bottom Wall Bottom Wall Bottom Wall Bottom Wall Bottom Wall Bottom Wall 


Figure 12: GFSSP Network Model of Driven Square Cavity 
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Figure 13: Shear Driven Square Cavity Centerline Velocity Distribution 

The predicted velocity field and pressure contours are shown in Figure 14. The recirculating 
flow pattern and stagnation of flow near the top right corner are clearly shown in the figure. The 
predicted stream traces from calculated velocity field is shown in Figure 15. 
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Figure 14, Predicted Velocity Field and Pressure Contours 
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Figure 15. Predicted Stream Traces in the Driven Cavity 

CONCLUSIONS 

The numerical algorithm of GFSSP has been extended to calculate multi-dimensional flow. 
GFSSP’s unstructured nodal network provides the ability to account for transport of scalar 
variables in n-dimensional space. Therefore, the formulation of the conservation equations for 
mass, energy and fluid species was not modified. However, the momentum equations are one- 
dimensional. This paper has described the extension of the momentum equations to calculate the 
shear force due to velocity gradient and transport of momentum due to transverse velocity. The 
implementation of these terms has been validated by comparing the numerical solution with 
three classical benchmark problems: Poiseulle Flow, Couette Flow and flow in a driven cavity. 
Excellent comparison has been obtained in all three cases. Future work will include turbulent 
flow and heat transfer. 
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NOMENCLATURE 


Symbol 


Description 


A Area (in 2 ) 

c Constant of Integration 
g Gravitational Acceleration (32.174 ft/sec ) 

m Mass Flow Rate (lb m /sec) 

m, n Number of Branches 
p Pressure (psi) 

t Time (sec) 

u Fluid Velocity in the x-Direction or in the Direction of Interest (ft/sec) 

v Fluid Velocity in the y-Direction or Normal to the Direction of Interest 

(ft/sec) 

V Volume (in 3 ) 

w Fluid Velocity in the z-Direction or Nonnal to the Direction of Interest 

(ft/sec) 

x, y, z The Cartesian Coordinates 
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Greek Symbols 

jli Fluid Viscosity (lb m / {ft sec}) 

p Mass Density (lb m /ft 3 ) 

0 Angle 

Subscripts 

i, j, k Summation Indices 
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